##########################################################
# John Henderson
# Replication Data for: "Hookworm Eradication as a Natural 
#   Experiment for Schooling and Voting in the American 
#   South", Forthcoming in Political Behavior, May 20, 2017
# 
##########################################################
#
#  figure_3.R
#  -- file produces qqplots in figure 3
#
##########################################################

rm(list=ls())

library(Matching)
library(AER)                           
library(foreign)

setwd('~/Dropbox/hookworm/replication')
source('funs/lower.bound.R')
            
load('data/hookwormData.Rdata') 

##########################################################  
# data required 
#  - mMat :: matching matrix used in matching and balancing
#  - oMat :: outcome matrix of county education and turnout rates
#  - zMat :: intervention matrix 
#  - pMat :: placebo matrix
# hookSouth :: combines all objects into one dataset, w/ additional covariates  
##########################################################  

##########################################################
##########################################################
# BEFORE MATCHING
##########################################################
##########################################################

d1=oMat$school1920
d0=oMat$school1910    

y0=rowMeans(cbind(oMat$vote1916,oMat$vote1920))
y1=rowMeans(cbind(oMat$vote1928,oMat$vote1932))

tr=zMat$z_campaign

d=d1-d0
y=y1-y0
d_dd=c(d0,d1)
y_dd=c(y0,y1)
tr_dd=c(tr,tr) 

xmat=mMat[,-c(1:3)][,1:19]
xmat=as.matrix(rbind(xmat,xmat))

year=c(rep(0,length(tr)),rep(1,length(tr)))
temps=as.data.frame(cbind(d_dd,y_dd,tr_dd,year))   

z_hook=zMat$z_hook_rate

{   
pdf(file="~/Dropbox/hookworm/replication/figures/figure_3a_before.pdf")
qqplot(z_hook[tr==1],z_hook[tr==0],ylab='Prior Hookworm Rate Without RSC Intervention',xlab='Prior Hookworm Rate With RSC Intervention')
abline(a=0,b=1,lty=2)    
ts=t.test(z_hook[tr==1],z_hook[tr==0])$p.      
text(cex=1.85,paste("p =", round(ts,digits=3),sep=' '),x=.2,y=.85)  
dev.off()
}
       

##########################################################
##########################################################
# AFTER MATCHING
##########################################################
##########################################################

mMat=bMat=mMat[,-c(1:3)][,1:19]

# pre-sets for 'lower.bound' fitfunc :: lower.bound only improves if total balance improves
floors=NULL					# :: allows to preset values at the floor overall
floor.values=NULL           # :: allows to preset values at the floor for each covariate
reweighted = -1 			# :: fits over a single summary measure of balance
exacts=array(0,ncol(mMat))  # :: exact matching for covariates 
exacts[1]=0    

starts=c(0.8639389,12.5338719,14.1059201,232.6227015,13.3530437,0.8680492,809.5765534,2.3393068,1.8030039,0.9325587,44.5785061,1.5006086,65.5576631,2.8099156,9.4050025,43.4617957,453.9672492,357.1535636,388.9175985)     


set.seed(1005)
genout=GenMatch(starting.values=starts,exact=exacts,MemoryMatrix=T,fit.func=lower.bound,Tr=tr,X=mMat,BalanceMatrix=bMat,M=1,pop.size=2,max.generations=5,wait.generations=3,hard.generation.limit=T,estimand="ATT",ties=T)

mout=(Match(Tr=tr,X=mMat,estimand="ATT",M=1,Weight.matrix=genout,exact=exacts,ties=T))

dtr=d[mout$index.treated]
dct=d[mout$index.control]
ytr=y[mout$index.treated]
yct=y[mout$index.control]
w=mout$weights
trt=tr[mout$index.treated]
tct=tr[mout$index.control]

ymat=c(ytr,yct)
dmat=c(dtr,dct)
zmat=c(trt,tct)
wmat=c(w,w)    

tt=tapply(trt,mout$index.treated,weighted.mean,weight=w)
cc=tapply(tct,mout$index.treated,weighted.mean,weight=w)
htt=tapply(z_hook[mout$index.treated],mout$index.treated,weighted.mean,weight=w)
hcc=tapply(z_hook[mout$index.control],mout$index.treated,weighted.mean,weight=w)

 
{
pdf(file="~/Dropbox/hookworm/replication/figures/figure_3b_after.pdf")
qqplot(htt,hcc,ylab='Prior Hookworm Rate Without RSC Intervention',xlab='Prior Hookworm Rate With RSC Intervention')
abline(a=0,b=1,lty=2)    
ts=t.test(htt,hcc)$p.value
text(cex=1.85,paste("p =", round(ts,digits=3),sep=' '),x=.2,y=.85)
dev.off()
}   

#END      